Bose-Einstein condensates on slightly asymmetric double-well potentials 
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An analytical insight into the symmetry breaking mechanisms underlying the transition from 
Josephson to self-trapping regimes in Bose-Einstein condensates is presented. We obtain expressions 
for the ground state properties of the system of a gas of attractive bosons modelized by a two site 
Bose-Hubbard hamiltonian with an external bias. Simple formulas are found relating the appearance 
of fragmentation in the condensate with the large quantum fluctuations of the population imbalance 
occurring in the transition from the Josephson to the self-trapped regime. 

PACS numbers: 



I. INTRODUCTION 

The dynamics of cold bosonic atoms in double- well po- 
tentials has deserved a great deal of attention in the last 
decade. See [l|, [2| for a careful review of the early find- 
ings, and [3[ for a more recent update. Of particular 
interest here are the seminal works of Smerzi et al. |J| 
and Milburn et al. '5(. The latter managed to derive a 
simplified many-body hamiltonian with semiclassical pre- 
dictions similar to those of |4|] but with the important ad- 
vantage of allowing the study of the quantum fluctuations 
on top of the semiclassical quantum averages. In both 
cases, the most remarkable feature reported was the ex- 
istence of a new phenomenon, macroscopic self-trapping, 
directly linked to the atom-atom interaction. Most of 
the studies focused on the case of repulsive atom-atom 
interactions. However, more recently interest has grown 
on the properties of these systems when this interaction 
is made attractive and in particular on the appearance of 
cat like ground states whose description goes beyond the 
usual semiclassical approximations [6|— lS|j . The structure 
of these ground states is determined by the ratio between 
interaction and tunneling strengths. 

For repulsive interactions, the Bose-Hubbard model 
has been widely used to study the transition between the 
Josephson and the self-trapped regimes 0, Is l9Hllj . It is 
also the natural choice for attractive interactions. Vary- 
ing the parameters of the model allows to easily scan the 
properties of the system. The transition from the Joseph- 
son to the self-trapped regime is directly related to the 
properties of the spectrum of the system as we increase 
the atom-atom interactions. In a recent manuscript the 
strongly correlated nature of the ground state of the sys- 
tem (for attractive interactions) in the transition region 
has been described [12|. This state is similar to the cat- 
like states described in Ref. [13, [lj| for the case of in- 
ternal Josephson-like behavior. Both studies find a deep 
relation between the appearance of a bifurcation in the 
semi-classical description and the existence of strongly 
correlated, cat-like, ground states in the systems. Sim- 
ilar states also appear in the nucleation of vortices in 
BECs [15| . The presence of the bifurcation points to the 



need of a fully quantum description that goes beyond the 
usual semiclassical approach. However the exact solution 
of the Bose-Hubbard model involves numerical calcula- 
tions that complicate the interpretation of the results. It 
is thus highly desirable to develop also simplified models 
that give analytical insights into the relevant physics. 

In this article we extend the previous work of Ja- 
vanainen et al. [6] , Jaaskelainen and Meystre [7j , and the 
more recent one of Shchesnovich and Trippenbach [8L Il6j| . 
Starting from Bose-Hubbard hamiltonians, these authors 
obtain approximate, Schrodinger-like, equations for the 
dynamics of the Fock space amplitudes of a system of N 
atoms on a double-well potential. We will first rederive 
these equations and, afterwards, will go one step fur- 
ther and obtain simple analytical solutions under certain 
premises. These are found to be in good agreement with 
the exact results and help to explain the delicate balances 
involved in the transition between the two regimes. 



II. QUANTUM MODELS FOR THE GROUND 
STATE 

The time dependent Schrodinger equation governing 
the evolution of the system of N atoms reads, 



-9 T |<i> >=#!$>, 



(i) 



with the (Bose-Hubbard) Hamiltonian written in reduced 
units as , 

H = -j^(a\a 2 + ala 1 ) + ^n 1 + ^{n 2 1 +n 2 2 ).(2) 

Where fii = a\ai, N is the total number of atoms in the 
system and ai (aj) is the annihilation (creation) operator 
for well i. The Fock state basis is written as |m, n^) and 



1 To ease the comparison with the work in Q we will make use of 
their notation. 



has N + 1 vectors, {|0, N),...,\N, 0)}. The two param- 
eters governing the dynamics are 7, which measures the 
ratio between the contact atom-atom interactions and 
the hopping strength, and e which is a symmetry break- 
ing bias, also divided by the hopping strength 2 . For this 
study we take e < which promotes |1). 

With this sign convention, 7 > and 7 < correspond 
to repulsive and attractive atom-atom contact interac- 
tions. 

The solution to © can be expanded in the Fock basis 
as, |$(i) >= 2fc=o c k{t) \k 7 N — k > , leading to an 
equation for the time evolution of the coefficients c k of 
the form d: 



i d 
l^~r c k = -(h-iCk-i + hck+i) + a k c k , (3) 



with b k = jjy/(k + l)(N-k), and a k = ^[k 2 + {N - 
k) 2 ] + jjk. Now, it is useful to introduce a new variable, 
X = k/N, and, assuming that h = 1/iV is small, we define 
\S?(x) — c k /VN and b(x) = b k - The next step is to make 
x a continuous variable 3 . This leads to Eq. (5) of [g], 

ihd T V(x) = - [e~ i§ b h {x) + b h (x)e ip ] tf (3) + a(x)*(x) , 

(4) 

where p = —ihd x , and 



b h (x) = [(x + h){\ - x)} 1 / 2 , 
a(x) = j[x 2 + (1-x) 2 } +ex. 



(5) 



Eq. (UJ) describes the time evolution of S&(x), which is 
a continuous interpolation of the Ck/yN- This interpo- 
lation assumes a certain degree of smoothness for the Cfc, 
so that for most k: sign(cfc + i) = sign(cfc). This is the 
case for the ground state with either repulsive or attrac- 
tive interactions, provided the symmetric state is lower 
in energy than the antisymmetric one (our case) 4 . 

From the formally exact Eq. (HJ, we first retain up to 
terms of order h 2 , 



e ±lp ~ 1 ± ha. 



1 



2o2 



h 2 d. 



(6) 



and 



b h (x) ~ b {x) + hd h b h (x)\ h=0 + (l/2)h 2 d 2 h b h (x)\ h=0 , (7) 



2 To compare to our previous work Il2l . lets note that, A = —27, 
and £R cf i 2 /J = e/2. 

3 Formally, we then have, bfc_i = e~ h9:c b(x), c^_i = 
T/Ne- hd *t>(x), and c fe+1 = v / A T e'» e *<I'(z). 

4 There are however other cases for which the cjj alternate in 
sign : CfcCfc + i < 0, Vfc. One example is the mean field state 
$ = [(l/ v / 2)(aj - al)] N \vac) As in @|, one way of dealing 
with this case is by introducing an auxiliary smooth function 
*(ik) = e mk V(x). 




FIG. 1: (color online) (a) The potential V for several values 
of 7 and zero bias (e = 0), (b) Separation between the two 
minima (solid-black) and height of the barrier as function of 
7, (dashed-red) , also in the case of zero bias. 



to get, 



ihd T i)(x) ~ ~h 2 b (x)ij"(x) - h%{x)ij)' {x) 

+ [o(a:) + V 2 (x)^(x) (8) 



where 



V 2 {x) = -2b {x) + h{b' (x) - 2d h b h (x)\ h = ) 

- h 2 (d 2 h b h (x)\ h=0 ~ d h b' h (x)\ h = + ~&o(z)) • 

(9) 

This equation is similar to Eq. (13) of Q if we retain 
only the h° term in T^(x). This approximation respects 
the symmetry and behaves well at the edges of the Fock 
space. Besides, taking into account that, 



P\/x(l - xjftip = -h 2 d x (b d x %lj) = -h 2 (b' n il>' + b tp") , 



(10) 



we arrive to 



ihd T tp = + [p^/x(l - x)p + a(x) + Vq(x)]i(j(x) , 



(11) 



with V (x) = -2 v /x(l -x). 

The range of x is the interval [0, 1]. It will be conve- 
nient to introduce the variable z = 1 — 2x, z G [—1,1]. 
In terms of the new variable, Eq. (fTTj) becomes 

ihd T tp(z) = ( - 2h 2 d z ^l - z 2 d z + V(zUi/>(z)(12) 

with V(z) = (l/2){ r yz 2 — ez) — yl — z 2 . This equation is 
a Schodinger-like equation with an effective z dependent 
mass: —h 2 /(2m) = — 2h 2 ^/l — z 2 . For convenience, we 
normalize ip(z) as 1 = /_. dz \i/j(z)\ 2 . 

As discussed in Refs. [16| and [13], the population im- 
balance is an appropriate order parameter. It is defined 
as, 

N 

I = ((*\n 2 - him/N = ^4 (1 - 2k/N) , (13) 

fe=0 
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FIG. 2: (color online) The potential V, black, compared to 
the three parabola used to approximate it near its extrema in 
building the analytical model. N = 50, 7 = —2.1. 



which in the continuous variable reads, / = 
j_ 1 dz\ip(z)\ 2 z , where we have used J2k=o c k = 1 an< ^ 
the normalization of the ^>{z). 

A. Effective two mode model 

We first study the solutions of Eq. ([T2"j) for the station- 
ary case with zero bias, and concentrate on the region 
where the classical bifurcation appears, 7 < — 1. 

For 7 < — 1, V(z) has two symmetric minima, which 
get deeper as 7 decreases, see panel (a) of Fig. [T] The 
separation between the two minima, 2^/1 — I/7 2 and the 
barrier height (A = (l + 7) 2 /(2|7|)) are depicted in panel 
(b). 

Since V(z) has pronounced minima and maximum we 
first determine the g.s. energy of the system by approxi- 
mating these minima by parabolas, see Fig. [3J 



V(z) ~ V{z m ) + -V"(z m )(z ± z m f 



(14) 



with z m = y/1 - I/7 2 , V(z m ) = 7/2 + 1/(27), and 
V"(z m ) = 7 — 7 3 . As an additional approximation, in 
the first term in the r.h.s. of Eq. (fT2]) we replace zby z m , 
since the wavefunction, ip(z), is sharply peaked around 
the minima at ±.z m . The stationary pseudo-Schrodinger 
equation can then be written as, 



-2^ 2 + i(7 4 



l 2 ){z-z m f 



i> 



= [-7^+(l + l/7 2 )/2]V- 



(15) 



and its spectrum can be readily obtained by comparing to 



the harmonic oscillator: h /(2m) 
1) and huj/2 



2/r, muj 2 -t 7 (7 



~^ l l\Jl 1 ~ 1- From this analogy one 
easily gets the ground state energy: 



E„ 



Yl + \ + ^ 



This approximation turns out to be very reasonable as 
it can be seen in Fig. [3] by comparing the g.s. energies 
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FIG. 3: (color online) Ground state energies of Eq. 1)120 . solid- 
black, compared to the approximate ones given by Eq. (|16l) . 
dashed-red, for N = 20, 50 and 90, respectively. 



obtained using Eq. dTrJl) and the ones obtained by solv- 
ing numerically the pseudo-Schrodinger Eq. ((T2|) . The 
approximation is particularly good for large N ( small 
h = 1/-/V), and 7 not too close to the bifurcation point, 
7 = — 1, as shown in the inset of the figure. 

For our analysis we assume that the harmonic oscilla- 
tor estimate of the ground state is accurate enough and 
that the shape of ip{z) around z m (and also — z m ) is that 
of the g.s. harmonic oscillator, i.e. gaussian with the 
parameters also determined from Eq. (1161) . These nor- 
malized solutions will be denoted \TV) or ^L^'°''(z), and 

\C) or "F^ ° (z). They correspond to the clockwise and 
counter-clockwise rotor solutions to the pendulum in the 
language of Ref. 4] . They form the "two-mode" basis for 
the model to be described below. Notice also that these 
two modes are built in the Fock space, they are useful 
to describe the system during the bifurcation and in the 
self-trapping regime. In fact, the model works soon after 
the bifurcation starts, once the overlap (C\1Z) is small 
enough. This overlap is given by, exp[2(7 2 — l) 3 / 2 /(7/1)]. 
As can be seen in panel (a) of Fig. HJ this overlap (black 
dot-dashed line) is essentially zero for 7 < —1.05 and 
N = 50. For these basis vectors, (£\z\C) — —z m and 
(1Z\z\1Z) = z m . In this two-mode model, an approximate 
solution can be written as: 



coscn|£) + s'ma\lZ) 



(17) 



In the absence of bias the ground state is the symmetric 
combination, and therefore a — it/ 4. The expectation 
value of z in the state (fT?]) is easily calculated, (z) = 
— z m (cos 2 a— sin a) — — z m cos(2a). The dispersion of z, 
a 2 = (tfj\z 2 \^) — (ip\z\ip) 2 , can also be readily computed. 
Retaining up to order h, 



2 2-2 

g, = z rn sin 



2a - /i/^vV - 1) • 



(18) 



(16) When no bias is present, 



= and a 2 = 



h/d^/l 2 — !)• The large values of a 2 are due to the 
fact that the ground state is cat-like, that is, the wave- 
function has two peaks in z. 



When a finite but small bias, |e| << 1, is considered, 
a is determined by the balance between the tunneling 
across the middle barrier of V and the bias term. In the 
absence of tunneling, \ip) will consist only of \C) or of 
1 72.) depending on the sign of the bias constant, e. In the 
latter situation the main c ontributor to a z is the other- 
wise small — /i/(7V7 2 — 1) term, as the sin 2 2a becomes 
negligible. 

To take into account both effects, the tunneling across 
the Fock space barrier and the bias term, we write an 
effective Hamiltonian for this two-mode model 



u = 



-t -ez m /2 



(19) 



The eigenvalues of W are E$ = —\Jt 2 + e 2 z 2 m /A = —Ea- 
Thus, one way to obtain the value of the tunneling con- 
stant, t > 0, is by computing the energy splitting be- 
tween the symmetric and anti-symmetric states in the 
double well in the absence of bias, t — (1/2)AEas = 
(1/2)(E A -E S ). 

To estimate this energy splitting in the double- well we 
will neglect the z dependence of the effective mass. This 
brings our problem to a Schrddinger-like equation and we 
can make use of the WKB approximation as developed 
in Razavy's book, (17l |. The energy splitting is written in 
Eq. (3.109) of that book: 



AE AS = — exp 

7T 



z + 



^J(2m/h 2 )(V(z)-E)dz 



(20) 



In our case we use the equivalences given above Eq. (|16[) 
to assign values to Ti /(2m) and hu>. The value of E 
corresponds to the ground-state energy and is taken 
from the harmonic oscillator approximation defined in 
Eq. (I16[) . To compute the tunneling integral, and obtain 
analytic results, we will approximate V(z) by an inverted 
parabola, V(z) ~ — 1+ ^(7+ l)z 2 , see Fig. [2] This gives 
an energy splitting, 



h/y 



AE AS = -2^V7 2 - 1 exp 

77 



■K \E\-\ 



2h ^/i^pi 



(21) 




FIG. 4: (color online) The three smaller plots above, (d), (e), 
and (f), show |V>( 2 )| 2 computed with the two mode model for 
7 = —1.6, —1.4, and —1.2 respectively, (a) Different ground 
state properties predicted by the analytic model for three dif- 
ferent regimes, separated by blue dotted lines. The three re- 
gions correspond to the bias-dominated fully-condensed (i), 
transition-asymmetric cat-like (ii), and symmetric cat-like 
(iii). We depict the behavior of — {z) g . s ., black solid line, and 
<j z , black dashed line, as a function of 7. z m is shown in thin- 
blue. The overlap {C\1Z) is the black dot-dashed line. The 
middle panel, (b), contains the condensed fractions, n± (up- 
per and lower black solid lines), and the mixing angle 0/2 of 
the ip+ eigenvector, see Eq. (J26J) , of the one-body density ma- 
trix associated to the ground state of the two-mode model. 
The lower panel, (c), depicts the energy splitting between 
the symmetric and antisymmetric (ground and first excited) 
states. The solid black line is calculated with Eq. (|22[1 . The 
blue empty circles are the exact numerical result obtained by 
solving Eq. (|12p . Red dot-dashed and red dashed lines depict 
2i and ez m respectively. In all cases e = —0.0002 and N — 50. 



This estimate of the splitting in absence of bias, is used 
to extract the value of t and it is in very good quantita- 
tive agreement with the exact energy splitting obtained 
by solving Eq. (fl~2"]) in absence of bias or as will be dis- 
cussed later, when the bias is not dominant. It is an 
important improvement on the estimate given in Ref. [8[ 
(cf. Eq. (31)). The value of 2i is plotted (red dot-dashed 
line) in panel (c) of Fig. 2] 

From the two mode Hamiltonian (Eq. ([T9"))) we find 
that the predicted "full" (including the bias) symmetric- 
antisymmetric energy splitting is 



AE 



' { 2s ) =2^/(ez m /2) 2 +t 2 



with two well-defined limits, 2i and sz m (red dot-dashed 
and red dashed lines of panel (c) of Fig. 2] respectively), 
depending on whether tunneling or bias is the dominant 
contribution. The ground state can be readily computed, 
Hips = Esips, with -0s = (cos a, sin a) and tana = £ + 
Ve + 1 , where, £ = (ez m /2t). 

The one body density matrix for the Bose-Hubbard 
Hamiltonian can be written as, 



P= N 



[a\ai) (a\a 2 ) 

;4«i) (4«2; 



(23) 



(22) with Trp = 1 and (/(*)) = fdz\^(z)\ 2 . In the large-N 



model we have, a\a\ = (l+z)/2,a\a,2 = (l—z)/2, a\a 2 = 
a 2 a\ = \/l — z 2 /2. In the two mode model retaining up 
to order h 1 we get, 



P = 



with eigenvalues, 



n± 



1 (l + z m cos 2a y/l 



V^ 



1 — z m cos 2a 



(24) 



l±Vl-^sin^(2a) 



i±v/T^ 



O(fc). 



(25) 



The last equality directly relates the appearance of frag- 
mentation and the quantum fluctuations of the popula- 
tion imbalance measured by er z , and holds also in the 
transition region. Denoting the corresponding eigenvec- 
tors by: (p± — (cos8±/2, sin6*±/2) one finds 



tan — = 
2 



z m cos 2a =F \ 1 



.sin 2 (2a) 



V^ 



(26) 



B. Results 

The quantitative predictions of the analytical model 
are presented in Fig.|H The model provides a simple and 
yet deep understanding of the problem. For 7 < —1.5, 
region (i), the ground state of the system is completely 
asymmetric (see the left upper smaller plot of Fig. |3J), 
with a large population imbalance, i.e. it is localized 
on the well promoted by the bias (|£) for e < 0). The 
location of the wells in the Fock space is measured by 
z m . The blue thin line in panel (a) shows z m , which is a 
decreasing function of 7 and reaches zero at the bifurca- 
tion point, 7 = —1. In this regime, the bias dominates 
completely the hamiltonian of Eq. (fTT)|) , which therefore 
is essentially diagonal with eigenvectors \C) and \R). In 
these cases, a z is small and given by the spread of the 
occupied mode. In this region, the system is fully con- 
densate. Correspondingly, as shown in the panel (b), the 
eigenvalues of the one-body density matrix of the ground 
state, are 1 and 0. In the panel (c), it is clearly shown 
that 2i is very small compared with ez m , and the split- 
ting E^g (Eq. [21]) in this region is given by ez m . 

In region (ii) both the bias term ez m and the tunneling 
matrix element, t, are of comparable size (see panel (c)). 
In that region the ground state is an asymmetric cat-like 
state (see the central upper smaller plot). The value of 
-(z) is decreasing and has large fluctuations, shown in 
the region (ii) of panel (a), mostly due to having both 
modes populated simultaneously. The system is frag- 
mented into two condensates, i.e the eigenvalues of the 
one-body density matrix are both different than zero, 
being ip + more populated, see panel (b) of the figure. 
The macro-occupied state, tp+, varies from close to |1) to 



(l/v / 2)(|l) + |2)) as can be seen from the behavior of the 
mixing angle + /2 (blue dashed line in panel (b)). 

The tunneling term, t, grows exponentially and is re- 
sponsible for the abrupt decrease of z as one crosses 
into the region, (iii). There, the two mode Hamiltonian 
is dominated by the tunneling matrix element, and the 
wave function becomes symmetric-cat-like, as shown in 
the right upper smaller plot. The cat like nature of the 
state reflects in a small value of the imbalance and a 
sizeable a z . Notice that when 7 approaches -1, the two- 
mode approximation is not valid anymore, because the 
(C I TV) , shown by the black-dot-dashed line in panel (a) 
is not longer close to zero. 

The simplifications of the two mode model capture to a 
very large extent the features of the initial large-N model, 
Eq. (TT2"|) . as can be seen from the comparison of the exact 
splitting (blue empty circles) between the ground and 
first excited state and the two- mode one (black solid line) , 
shown in the last panel of the figure. 



III. CONCLUSIONS 

For attractive interactions, we have discussed here the 
appearance and properties of the cat states in the tran- 
sition from the tunneling to the interaction dominated 
regimes. The semiclassical approach shows that at some 
critical value of the ratio of atom-atom to tunneling en- 
ergies, 7, a bifurcation appears. We have constructed an 
effective "two mode model" built from the solutions cor- 
responding to the branches of the bifurcation. And have 
shown that it is quite successful in describing not only 
the average imbalances in the presence of bias but also 
the quantum fluctuations. The model includes the main 
quantal effects since its predictions are in good agreement 
with exact solutions of the Bose Hubbard hamiltonian. 

Our starting point is an already developed large-N 
model [6J, |8| which gives accurate numerical predictions 
for N > 30. We have managed to obtain analytic ex- 
pressions valid on the transition region from self-trapped 
(bias-dominated) to the symmetric cat-like region (dom- 
inated by the tunneling term). The main feature is the 
sharp change of the population imbalance as 7 is de- 
creased. It is due to the delicate interplay between the 
bias and an effective tunneling term with a sharp depen- 
dence on 7. 

The fact that Eq. (|12p captures most of the non-linear 
dynamics of the original Bose-Hubbard hamiltonian im- 
plies that the complex dynamics of systems governed by 
non-linear terms can be studied very conveniently with 
pseudo-Schrodinger equations governing the dynamics in 
Fock space. In these equations the main effects of non- 
linearity appear through a barrier-like potential. In this 
way, some of the dramatic consequences of the non- 
linearities are mapped into exponentially varying mag- 
nitudes in linear systems. We have exploited this advan- 
tage to introduce further approximations and construct 
a very simple effective two-mode model to reproduce and 



explain the main properties of the transition. 
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